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rvq We analyse a simple 'Stokesian squirmer' model for the enhanced mixing due to swim- 

ming micro-organisms. The model is based on a calculation of Thiffeault & Childress 

' -h' [Physics Letters A 374, 3487 (2010)], where fluid particle displacements due to inviscid 

^..^ swimmers are added to produce an effective diffusivity. Here we show that for the viscous 

■■^ case the swimmers cannot be assumed to swim an infinite distance, even though their 

total mass displacement is finite. Instead, the largest contributions to particle displace- 
ment, and hence to mixing, arise from random changes of direction of swimming and 
are dominated by the far-field stresslet term in our simple model. We validate the re- 
. N suits by numerical simulation. We also calculate nonzero Reynolds number corrections to 

• i-H the effective diffusivity. Finally, we show that displacements due to randomly-swimming 
^^ squirmers exhibit PDFs with exponential tails and a short-time superdiffusive regime, as 
r^ found previously by several authors. In our case, the exponential tails are due to 'sticking' 
Oh near the stagnation points on the squirmer's surface. 

(N 
> 

O 

■^ 1. Introduction 

. Swimming creatures affect their environment in many ways, and one which has received 

attention recently is how they mix the surrounding fluid. This phenomenon is called 

biogenic mixing or biomixing. The most striking and controversial setting is the ocean: 

^^ Dewar et al. (2006) suggested that marine life might have an impact on vertical mixing 

^_H in the ocean. 

ILJ Katija & Dabiri (2009) proposed that the dominant effect involved in biomixing is the 

• ^H mass displacement due to a swimming body. This phenomenon is called Darwinian drift, 
^ after Darwin (1953), though the displacement due to a moving cylinder was obtained 

by Maxwell (1869). Thiffeault & Childress (2010) derived the effective diffusivity of an 
'ideal gas' of randomly-distributed non-interacting swimmers, and showed that it depends 
on the induced squared-displacemcnt of fluid particles by the swimmers, as opposed to 
the net mass displaced. 

In the present paper we apply the techniques of Thiffeault & Childress (2010) to 
swimmers in the low- Reynolds number regime. Since the experiments of Wu & Libchaber 
(2000), there has been considerable work in that direction. Two general types of swimmers 
arise: pushers or pullers, depending on whether the propulsion takes place ahead or behind 
the swimmer's centre of mass. Several authors (Dombrowski et al. 2004; Hernandez-Ortiz 
et al. 2006; Saintillian & Shelley 2007; Underbill et al. 2008) have pointed out that at 
high volume fractions the pushers align spontaneously with their neighbours, producing 
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Figure 1. Definition of impact parameters a and h, displacement A = /S.\[a, b), and swimming 
patli lengtli A. In tliis picture tlie parameter b is positive; negative b corresponds to tlie swimmer 
starting its trajectory past tlie point of initial closest approach. The filled dot is the initial 
position of the target particle and the hollow dot is its final position after the swimmer has 
moved by a distance A. The 'interaction disk' of radius R is also shown. 

large-scale motions. At low volume fractions, pushers and pullers behave similarly, mov- 
ing in a relatively uncorrelated manner. It is this limit that we treat in this paper. We 
use Lighthill (1952) and Blake (1971) 's 'squirmer' model, where a swimmer is modelled 
as a sphere in Stokes flow with a prescribed tangential velocity distribution. These have 
recently been studied by Ishikawa & Pedley (2007) and Drescher et al. (2009) as models 
of the green alga Volvox. Squirmers could also be appropriate for Chlamydomonas (Lep- 
tos et al. 2009; Guasto et al. 2010; Drescher et al. 2010), as the observed velocity fields 
resemble those of the models described here. We find a formula for the effective diffusion 
coefficient due to the squirmers, and show that it depends predominantly on the far-field 
displacements imposed by the swimmer. Moreover, the dominant contributions arise from 
the turning motion of the swimmers, in sharp contrast to the potential fiow case (Thif- 
feault & Childress 2010) where the path length of straight swimming could be considered 
infinite for all practical purposes. We also show that the probability distribution function 
of displacements due to squirmers has exponential tails, as in Leptos et al. (2009), and 
that it exhibits a short-time supcrdiffusive regime as in Wu & Libchaber (2000). 

The formula we derive for the effective diffusivity has a wide applicability beyond 
squirmers, for instance to swimmers at larger Reynolds numbers with wakes. In that 
spirit, we derive finite (but small) Reynolds number corrections to the effective diffusivity 
using the Oseen form for the far field, and find that the effective diffusivity falls off as 
^g-0.61 f].Qjjj j^g Stokes (zero Re) value. 



2. Computing the effective diffusivity 

We first consider a single swimmer in two or three dimensions that moves in a straight 
line at constant speed U for a distance A. For simplicity, the swimmer is assumed axially 
symmetric with characteristic size £, and moves along its axis of symmetry. We are 
concerned with the effect of the swimmer on a target Lagrangian particle. There are 
two 'impact parameters' necessary to describe such a situation: the initial perpendicular 
distance to the particle, a ^ 0, and the relative distance b of the start of the trajectory to 
the point of initial closest approach (see figure 1). Note that b may be positive or negative, 
and is positive if closest approach would occur at a positive time. After the swimmer 
traverses its full path length A, the target particle is displaced by a distance Ax{a, b). 



Stirring by squirmers 3 

Inspired by Einstein's theory of Brownian motion (Einstein 1956) and following Thif- 
feault & Childress (2010), we now model the displacement of the target particle af- 
ter M encounters with swimmers as a series of 'kicks', each of which moves the particle 
by A\{ak, bk) in some random direction r^. The position xm of the particle is 

M 

Xm ^XQ + ^l::^x{ak.,bk)rk (2.1) 

fe=i 

and we will assume Xq = 0. The swimmers are assumed identical and non-interacting, and 
the encounters are taken to be statistically independent. Their directions of swimming 
are assumed isotropically distributed. In this case, the mean square displacement after 
M such encounters is given by 

M 

{\XM?) = ^(A^(a,,6,)) = M(A^(a,6)), (2.2) 

because the a^ and bk are identically distributed for each encounter, with the brackets 
denoting an average over a and b. (We could also average over distributions of A and C/, 
but we assume constant values for simplicity.) Our theory will apply in the dilute limit, 
which we will define more carefully below. 

The displacement A^ is assumed small compared to A, so that we can regard the target 
particle as approximately fixed during an encounter. We assume that each swimmer 
advances by a fixed length A and then changes direction randomly, so after a time t, 
Ut/X such direction changes have occurred. For the 2D case, the number of swimmers 
in a ring of radius L and thickness dL around the target particle is 2-KnLdL, where n 
is the number density of swimmers. To define an encounter, we introduce an 'interaction 
disk' (sphere in 3D) of radius R around the target particle (see figure 1), and consider an 
encounter as occurring when a swimmer enters that disk. If a swimmer does not enter the 
interaction disk, then the motion of the target particle is negligible. The interaction disk 
allows us to be more precise about what we mean by a dilute suspension of swimmers: 
we require nE? <C 1 {nR^ <C 1 in 3D) so that encounters are isolated from each other. 

The expected number M of such encounters after a time t is best written in terms of 
integrals over a = |isin0| and b = LcosO: 

2nUt /-^ /-A+ys^^^ ^2 



M :^ ^-- / / dbda = 2nRUt + TmUt^ . (2.3) 

For large path length A 3> i? the final term can be neglected, and we obtain the same 
expression as in Thiffeault & Childress (2010), where an infinite path length was assumed: 
for large path length, the frequency of encounters is independent of A. 

Next we compute the average squared displacement (2.2) by integrating A^- 



A- 



2nUt 



R /-A+V-R^-a^ 



M{/^i) = —- / Al(a,6)d6da. (2.4) 

^ Jo J-VW^^ 

We assume that /S.\{a, b) decays rapidly for large a^ £, where £ is the characteristic size 
of a swimmer. In that case we can replace i? ^ ^ by infinity; after also replacing Xm 
by x{t) we obtain 



TT+ pQO pQO 

{\x{t)\^) = —^ / Al{a,b)dbda^:4Kt (2D), (2.5) 

A Jo J-oo 

where we introduced the effective diffusivity k. In 3D the analogous expression can easily 
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shown to be 

{\x{t)\^) = ^^^ r r a^l{a,b)Abda=:QKt (3D). (2.6) 

Note that the 'interaction disk' of radius R has disappeared from (2.5) and (2.6): be- 
cause of the rapid decay of l^\{a^h) with a, we are free to overcount encounters when 
computing (|a;(t)p) because faraway encounters hardly displace the target particle. 

Thiffeault & Childress (2010) discussed models based upon potential flow past a cylin- 
der or a sphere and made use of an approximate form of (2.5)-(2.6). The approximate 
form effectively assumes that the particle doesn't move unless the swimmer actually 
crosses the point of initial closest approach (a good approximation for potential flow, see 
figure 4(a)). Thus /S.\{a^b) takes the 'top-hat' form 

A.(a,i) = |„^<°'- ■"<'<^' (2.7, 

I 0, otherwise. 

This approximation works very well for potential flow. With the approximate form (2.7), 
the effective diffusivities defined by (2.5)-(2.6) become 

{\Unj^^\a)da, (2D), 

\ft/n/„°°aA2(a)da, (3D), ^ ' 

as in Thiffeault & Childress (2010). Because the potential flow approximation is inappro- 
priate for micro-organisms, we will not use (2.8) but rather take the full form (2.5)-(2.6). 
Note that the number density scahng of (2.5)-(2.6) was derived by Underbill et al. 
(2008) using the Green-Kubo formula, but here we get explicit values for the prefactor. 
As we will see below, this allows us to determine where the dominant contribution to k 
arises from. In addition, the dimensionless prefactor can vary significantly depending on 
the swimmer under consideration. 



3. Stirring by squirmers 

Squirmers are a simple model of swimming micro-organisms in the Stokesian regime. 
Lighthill (1952) first introduced the model and Blake (1971) later extended the analysis. 
Recently, the squirmer model has been applied to describe the behaviour of slow-moving 
ciliates (see for example Ishikawa et al. 2006; Ishikawa & Pedley 2007; Drescher et al. 
2009), whose motion is propelled by arrays of pulsating cilia on their surface. A squirmer 
is generally a three-dimensional sphere moving in a Stokesian fluid with velocity field 
prescribed on its boundary. The velocity field generated by such a squirmer is derived 
by an expansion of solutions obtained by separation of variables. The steady swimming 
imposes a force-free condition between the fiuid and the body which in turn translates 
to algebraic equations for the expansion coefficients of the solution. 

Following Ishikawa et al. (2006), we choose a specific instance of an axially-symmetric 
squirmer of radius I, swimming along the positive z-axis at a constant speed U . Its 
free-space, steady, axisymmetric streamfunction in a comoving reference frame is 



^(,,.) = ip^L/ -i + ^ + ^^£:f ^-1 (3.1) 



in cylindrical coordinates where r — \J p^ + z^ — \J x^ + y"^ + z"^. The stresslet coefficient 
^ is a free, dimensionless parameter and the special case /3 = corresponds to a sphere 
in potential fiow; we shall take /3 = 5 in numerical examples as in Ishikawa et al. (2006). 
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Figure 2. Streamfunction contours for a steady squirmer with /? = 5, in a comoving reference 
frame. Note the wake or closed 'bubble' behind the squirmer. 





b/X = l 

Figure 3. Typical trajectories of a target particle in the fixed lab frame for different values of 
a and b due to a single squirmer moving from left to right, with I — U = 1, X = 100 and j3 = 5. 
The trajectories are offset vertically for clarity. The initial positions of the particle are marked 
by solid dots and the final positions by hollow dots. In each panel, solid dots from top to bottom 
correspond to \og{a/£) — 2, 1, 0, —1, —2, —3. 



Figure 2 shows the contours of the streamfunction (3.1) for a squirmer swimming from left 
to right, and figure 3 shows typical trajectories of a passive target particle displaced by the 
squirmer with different values of the impact parameters a, b. These types of trajectories 
were discussed by Dunkel et al. (2010) in the context of swimming organisms. 

Compared with the fore-aft symmetry of the potential flow case (3 — 0, a. squirmer 
with /3 > has a downstream wake or 'bubble', which we will discuss in section 4. In 
addition, a passive particle displaced by the squirming motion experiences a much larger 
excursion when the swimmer is still more than several body lengths away from the point of 
initial closest approach, while in potential flow this contribution is negligible to the total 
displacement. This suggests that the finite-path displacement function A>,(a, b) due to a 
squirmer has a more complicated dependence on the details of the swimming path through 
the values of a, 6, and A. A direct consequence of this feature is that the diffusivity integral 
can no longer be reduced to one dimension by restricting the integration domain of b as 
in equation (2.7). Instead, the full two-dimensional integral (2.6) has to be computed. 

Thiffeault & Childress (2010) showed that the integrals (2.5) and (2.6) for potential 
swimmers are dominated by "head-on" collisions, namely, significant particle displace- 
ment only occurs when the impact parameter a < ^ and when the swimmer passes 
by the particle. We next investigate the dominant contribution of the integral for the 
squirmers (2.6). Transforming the integral as 



x-H- 



aAl{a,b)dadb 



rVAA(a, 6) dlog(a/^) d(6/A), 



(3.2) 



we compare in figure 4 the spatial distributions of the dimensionless integrand i~'^a'^A\ (a, b) 
due to a sphere in potential flow and a squirmer in Stokes flow, for £ = U = I and A = 100, 
in the \og{a/£)-{b/X) plane. The far-field stresslet-only approximation to the displace- 
ment for a squirmer is also plotted as an inset in figure 4(b). The inset plot is almost 
identical to the full plot for any A » £. Note that there is a vertical strip near 6 = in 
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Figure 4. Integrand £~^o^A^(o,6) from (3.2) due to (a) potential flow past a sphere; (b) a 
squirmer with /? = 5; the near-field potential-flow pattern can be dimly seen for log a < 0, but 
is much weaker than the far field contribution. Inset: integrand corresponding to the far-field 
approximation of the squirmer flow. In all figures, i. = U = 1 and A = 100. 



6/A = 





b/\ = 0.5 



b/X = 1 



Figure 5. Trajectories of target particle for a/(. = 1.1 and different values of 6. The filled dot is 
the initial position of the target particle and the hollow dot is its final position after the swimmer 
has moved by a distance A = 100 to the right. The net displacement is small at 6 = A/2 because 
of cancellation as the particle moves towards and then follows the swimmer. 



each figure which corresponds to the situation where the initial particle position is inside 
the swimmer or its wake and is therefore excluded in the calculations. We shall return 
to the effect of particles trapped in the wake in section 4. 

For a sphere in potential flow, the dominant contributions to (3.2) come from a rect- 
angle (figure 4(a)) where the values of the integrand arc nearly uniform in h for h £ [0, A]. 
This justifies the approximation (2.7). However, for a squirmer the dominant contribu- 
tions to the integral (3.2) come from a much larger and irregular region with a/£ > 1. 
The contribution from this region is negligible in the potential case. 

What explains the strong peaks in figure 4(b) along the lines b/X = 0, h/X — 1? Figure 5 
provides the answer: a target particle with h/X = 1/2 undergoes a sizable excursion, but 
has a small net displacement since it moves towards and then follows the squirmer, giving 
back a large part of the displacement it first underwent. In contrast, for b/X — Q and 
b/X = 1 the squirmer starts or stops (i.e., changes direction) before the particle reverses 
direction and gives back its displacement. 

Now that we know where the contributions to (3.2) arise from, we apply the basic 
theory to the squirmer model and find 



K = 57.0[/nr, A>£ 



(3.3) 



which is more than 200 times bigger than the effective diffusivity due to potential swim- 
mers under the same set of parameters. One source of this large enhancement comes 
from the free parameter /3, equal to 5 here. Since at far field the dominant perturbation 
to the uniform flow in the stream function (3.1) is multiplied by /3, we expect that for 
large /9, k ~ 2.ip^Un£'^ (the far-field value, with coefficient fitted in figure 6(a)), whereas 
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Figure 6. (a) Dependence of efTective diffusivity on /? for A 3> ^: the solid line is the numerical 
integration of (2.6) with Aa computed with the streamfunction (3.1). For large A, k is inde- 
pendent of A. The dashed line is the effective diffusivity due only to the stresslet term, (b) The 
mean-squared displacement (solid line) of a target particle for 2 x 10® realisations of 10 squirm- 
ers, with £ = U = 1. The dashed line is the squared-displacement predicted by (3.3), using a 
number density n — 10~*£~^. 

the difFusivity converges to the potential value 0.266 Un£'^ as /? — >■ 0. Figure 6(a) shows 
how different values of /? change the effective diffusivity k. Using the far-field form of our 
model, we can recover the diffusivity 1.53 /im'^/sec measured by Leptos et al. (2009), for 
volume fraction 1.6 %, velocity U = 100 /xm/sec, and size (, = b /im, using a stresslet 
amplitude (3 = 0.6. The linearity of the diffusivity in the velocity and volume fraction 
(for dilute suspensions) has been pointed out by Underbill et al. (2008) and Leptos et al. 
(2009), but is verified here from first principles with computed prefactors. 

Next we validate the results with direct simulations conducted as in Thiffeault & Chil- 
dress (2010). For each trial, we use 10 randomly-distributed, non- interacting squirmers 
with i — U = 1 moving in random initial directions in a large periodic box of size 1000. 
Each squirmer moves in a straight line for a length A — 100 and then turns in a new 
random direction. The turning times are randomized so that the squirmers do not all 
turn simultaneously. The displacement statistics are identical when this asynchrony is 
relaxed. The target particle initially at the origin is displaced by the superposition of the 
flows created by each squirmer (this is a good approximation at such low densities) . Fig- 
ure 6(b) shows the mean-squared displacement (|a;p) of a target particle over 2 x 10® trials 
(realisations). The agreement between the solid curve (numerically-simulated (|a;p)) and 
the dashed line (6Kt with k, computed with the integral (2.6)) verifies the the theoretical 
prediction. We also simulated a large number of particles advected by a single realisa- 
tion of 30 squirmers: the squared-displacement averaged over all the particles obeys the 
theory for large times, as expected. 



4. Wake transport 

In the previous section we have ignored the wake or closed 'bubble' behind the squirmer, 
visible in figure 2, since it is helpful to treat it separately. Assume that a particle is 
initially trapped in the swimmer's wake: it will then be displaced by the full swimmer's 
path length A, until the swimmer changes direction and new particles are introduced in 
the wake. We can thus approximate the displacement Aa to be equal to A if the particle 
is inside the wake. The wake contribution to the effective diffusivity is then obtained 
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Figure 7. The role of the Reynolds number Re. (a) Integrand l'^a^A\{a,h) from (3.2) for 
displacements due to the far-field Oseen flow (5.1), for Re — 0.5 and A — 100. (b) The effective 
diffusivity k from (2.6) as a function of Re. For small Re we recover (3.3). 



from (2.6) as 



TT r 

'«wako = 'TTT^^ I 27ra dadb = \UnV^ 

"^ J wake 



wake 



A 



(4.1) 



where Vwake is the volume of the wake 'bubble'. Note that nV^wake is the total volume 
fraction of wake bubbles in the fluid. The wake effective diffusivity (4.1) is potentially 
much larger than that due to displacements outside the wake bubble, since (4.1) is directly 
proportional to A. In practice, many swimmers do not exhibit a wake bubble, and even 
for those that do the time-dependence of the velocity field near the swimmer and the 
effect of other swimmers can cause particles to jump in and out of the wake, which may 
lower (4.1) considerably. Molecular diffusivity will also limit the path length over which 

material remains trapped in the wake: for A > ^^ake/'^moiccuian significant amounts of 
material will have leaked out. 



5. Nonzero Reynolds number 

In the previous section, the basic theory was applied to model squirmer-induced mixing 
in Stokes flow, for which the Reynolds number Re = and inertial effects are neglected. 
We identified the dominant contributions as arising from the far field (stresslet). To 
generalise the squirmer far field and to account for finite (but small) Reynolds number 
effects on diffusion, we consider the analogous problem for the Oseen equations, which 
linearise the inertial forces (Oseen 1910). The singularity analogous to the squirmer far- 
field has streamfunction 



r^ \ Re r \ 



-Re (r-z)/2£ _ 1 



(5.1) 



It is easily checked that a,s Re —> this recovers the far-field displacement due to a 
Stokesian squirmer. Using this streamfunction we obtain numerically A\{a, b) for different 
values of Re and summarise the results for the effective diffusivity in figure 7(b). For 
small Re, we recover the Stokes value of the effective diffusivity; for larger values of Re, 
the effective diffusivity drops off as Re~°'^^ due to smaller particle displacements. 
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Figure 8. (a) Growth of mean-squared displacement {(Aa;)^) with time, with superdifFusive 
(short-time) and difTusive (long-time) regimes, (b) Rescaled a;-displacement distributions in the 
two regimes, showing the collapse after scaling. Inset: unsealed PDFs with exponential tails (6.1). 



6. PDF of displacements 

We close with some comments regarding the typical size of displacements imparted by 
the squirmers. In a recent letter, Leptos et al. (2009) measured experimentally the prob- 
ability distribution function (PDF) of displacements of a target particle due to repeated 
'kicks' by swimming micro-organisms. They observe that, unlike in Einstein's theory of 
Brownian motion, the PDF is non-Gaussian with clear exponential tails. We now show 
that squirmers also exhibit such tails, due to the stagnation points at their surface. 

The large displacements for squirmers are due to particles that lie very close to the 
axis of motion, i.e., with a <S^ £. For small a, the displacement function can be well- 
approximated with a 'top-hat' function as in (2.7). The probability of A\{a,b) being 
larger than a given value S is then 



P(Aa >S) = 



1 



V 



(da. 



Tradadb — - 

'Ax{a,b)>S ^ JA{a)>S 

For a <ti £, the squirmer displacement A (a) diverges as C£\og{i/a), where C is a numer- 
ical constant. Hence, the contribution for large S comes from small a, and we have 

A /■ . ttA 



P(Aa > 5) 



V 



(6.1) 



which exhibits exponential tails, as in Leptos et al. (2009). A simple argument using 
large-deviation theory shows that the PDF of a superposition of such exponential tails 
retains the same exponential tails. 

To make contact with Leptos et al. (2009), we statistically analyse displacements for 
our biomixing model via particle simulations. We examine the displacement histograms 
for different sampling intervals At up to 90, over which each swimmer travels a distance 
UAt. Figure 8(a) shows the mean-squared displacement ((Aa:)^) as a function of time: 
for short times, this grows as (Ai)^'*, in good accordance with the supcrdiffusivc regime 
observed by Wu & Libchaber (2000): they found an exponent between 1.5 and 2. For 
longer times, we recover a linear diffusive regime, as expected, with slope given by our 
theoretical prediction 2k. Figure 8(b) shows that the incremental displacement probabil- 
ity distribution along the cc-axis, P(Ax, At), rescales to a single distribution function in 
terms of the similarity variable Ax/^y{{Ax)^) in each of the two scaling regimes. How- 
ever, note that the exponential tails (6.1) are independent of At, and are only visible 
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in simulations after integrating for a long time (see inset in figure 8(b)). Thus, it is not 
clear if these are related the the exponential tails observed by Leptos et al. (2009) . 

The authors are grateful to Keith Moffatt for suggesting the wake transport calculation, 
to Eric Kunze for pointing out that molecular diffusion would limit the wake transport, 
and to the Institute for Mathematics and its Applications (supported by NSF) for its 
hospitality. SC was supported by NSF under grant DMS-0507615, J-LT under grant 
DMS-0806821. 
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